Resolving and correcting for kinetic biases on methane seep paleotemperature using carbonate ∆47/∆48 analysis

Methane-derived authigenic carbonate often constitutes the sole remaining record of relic methane seeps. The clumped (∆47) and oxygen isotopic composition of seep carbonates often yield inaccurate temperatures, attributed to kinetic isotope effects and modification of seawater isotope composition by hydrate water. Here, we analyzed the dual-clumped isotope (∆47/∆48) composition of authigenic carbonate from a modern methane seep. We demonstrate that aragonite forms closest to isotopic equilibrium such that its ∆47 can directly yield the correct formational temperature, whereas calcite is unambiguously biased by kinetic isotope effects. Numerical models show that the observed bias in the isotopic composition arises from rate-limiting dehydration/dehydroxylation of HCO3− alongside diffusive fractionation, which can be corrected for with analysis of carbonate ∆47/∆48 values. We demonstrate the utility of dual-clumped isotope analysis for studying seep carbonates, as it reveals the origin and magnitude of kinetic biases and can be used to reconstruct paleotemperature and seawater δ18O.

The PDF file includes:

Supplementary Text Figs. S1 to S9 Table S1 Legends for data S1 and S2 References
Other Supplementary Material for this manuscript includes the following:

Supplementary Text Setup and governing equations in the COAD-MS methane seep box model:
The COAD-MS model and the scripts used to generate Figures 2 and 3 can be accessed via the following Zenodo DOI: https://doi.org/10.5281/zenodo.11066480 The COAD-MS model simulates the chemical reactions governing dissolved inorganic carbon, oxygen, calcium, methane, sulfate, sulfide, total alkalinity and the carbonate minerals precipitating from the fluid.The default initial state of the model (except where otherwise specified) is an approximation of bottom water in the South China Sea, shown in Table S1.

Reaction equations and rates for methane -molecular oxygen -sulfate -sulfide
To model the sulfate-driven anaerobic oxidation of methane (SD-AOM) and aerobic oxidation of methane (AeOM) and the oxidation of sulfide, the following reactions (also shown in main text) are modeled in addition to those used in the original Watkins and Devriendt (2022) COAD model (26).
To summarize, for one mol of methane consumed, SD-AOM yields 1 mol of DIC and 2 mols of alkalinity and AeOM yields 1 mol of DIC an no alkalinity.The aerobic oxidation of sulfide consumes two mols of alkalinity, thus SD-AOM and subsequent oxidation of sulfide is essentially the same as simply aerobically oxidizing the methane from the perspective of carbonate saturation state and pH.The precipitation of CaCO3 consumes two mols of alkalinity and 1 mol of DIC.The effects of these processes on solution pH and CaCO3 saturation state are summarized in Figure S1.
The above chemical reactions are assumed to release CO2 and HCO3 -that are at isotopic equilibrium with respect to carbon, oxygen and clumped isotopes, defined by the fractionation factors of Beck et al. (59), Mook(60), and Hill et al. (61).This assumption is taken in order to test if there is a fundamental necessity for microbial AOM-specific fractionation to produce the measured isotopic disequilibrium of precipitating carbonates.There are arguments that these fractionation effects may exist, as it is understood that microbial processes can produce DIC that is out of isotopic equilibrium (46,62).We argue that there is value in initially assuming that no such fractionation occurs and testing this system purely from the perspective of DIC reaction kinetics responding to a chemical perturbation.Our findings show that the observed variance and co-variance in the δ 18 O-∆47-∆48 system can be readily described with this approach, demonstrating that there is no need for such additional reaction-specific fractionation factors.
A function, titled Seep_ODE.mhas been constructed in Matlab, which computes the instantaneous rate of change of CO2 (mass 44-48) Exchangeable Inorganic Carbon (EIC: HCO3 - +CO3 2-mass 60-64) as governed by the COAD model, as well as alkalinity, sulfate, sulfide and oxygen.The model uses the total concentration of dissolved inorganic carbon, and the alkalinity, to calculate the instantaneous pH using the CO2SYS model (57).Chemical transport with the boundary condition (seawater) is governed by the exchange rate, l, which is constant for all chemical components and is governed by Fick's law for diffusive transport, where Cexternal is either seawater or the methane supply.
Mass-dependent fractionation of isotopologues by this transport flux is allowed by varying l following a mass-dependent relationship, shown by example for a mass-63 isotopologue of dissolved carbon.
In a purely kinetic system, average particle velocities (and thus their rate of dispersal) at a given temperature are dictated by the square root of their mass, thus b would be 0.5; due to interactions between particles, such as the hydration of ions by surrounding water molecules, this number is typically much less.According to molecular dynamics simulations systems, b for DIC constituents can be up to 0.17, (28,45).In the models shown in the main text, a range of values is shown in Figure 3g-h between 0 and 0.15, however the default value in these models is 0.08, which appears to yield agreeable results with measured values.
As l increases, more reaction chemistry is able to occur before products and reactants are exchanged with the boundary conditions, resulting in variable extents of kinetic offsets.Models show four distinct stages during the progression.1) O2 consumption, with corresponding increase in ∑DIC and decrease in pH. 2) SO4 2-consumption, with corresponding increases in ∑DIC and pH, and offsets in isotope compositions.This occurs until 3) CaCO3 precipitation occurs, which is associated with a more limited increase in pH, as well as consumption of ∑DIC, Ca 2+ .4) Attainment of a steady state, where no further change in chemistry occurs and all fluxes are balanced.
When the output steady state values for carbonate flux, pH and precipitating carbonate are plotted relative to l, Figure S3 is generated.Figure S3 shows that there is an optimum rate of transport that allows carbonate precipitation to occur at an optimum rate, and that this optimum carbonate precipitation flux, also roughly coincides with the maximum kinetic offset in δ 18 O, ∆47 and ∆48.Crossplotting these output isotopic values produces the diagrams shown in Figure 3c Salinity affects a number of model behaviors, critically it is related to the saturation state of carbonate, as well as the rate of hydration and hydroxylation reactions.Increasing salinity results in lower equivalent saturation state of calcium carbonate, and thus higher concentrations of DIC are required to precipitate CaCO3, resulting in diffusive fractionation manifesting itself more in these systems.Figure S4 tests the COAD-MS model using salinities of 20, 30 and 40 ppt.Lower salinities than 20 ppt result in super-saturation of CaCO3 at the starting condition.
Temperature governs the equilibrium isotopic composition for δ 18 O, ∆47 and ∆48 of carbonate the model, however the rate constants for hydration/hydroxylation are both also related to temperature (64), thus warmer temperatures are associated with more rapid equilibration of the DIC pool, and thus with comparatively smaller kinetic offsets and different equilibrium states, this is shown in Figure S5, which tests the model at 0°C, 10°C and 20°C.The coefficients used for methanotrophy are not temperature-dependent, however if these are also increasing, then kinetic offsets may still be expected at warmer temperatures due to faster rates of SD-AOM and AeOM.
Similarly, the rate of equilibration and the saturation state of carbonate are both affected by solution pH, which can also vary in the modern ocean between 7.6 and 8.2, but may have ranged between 6.5 to 9 throughout geologic history (65).We explore this in Figure S6 with simulations run with pH values of 7, 7.5 and 8.It is unclear how realistic independently varying pH is, as it is generally coupled in seawater with Ca and DIC concentrations via the solubility product of calcite or CO2 removal and addition via photosynthesis or atmospheric exchange.Nevertheless, lower pH is associated with more rapid equilibration rates due to increased [CO2] and [H + ] concentrations facilitating faster hydration reactions, and thus lower pH is associated with less kinetic offsets than higher pH.Solution pH has a secondary effect related to carbonate saturation state, in that DIC concentrations must be higher at lower pH in order to result in mineral formation, a consequence of this is an enhanced contribution of mass-dependent fractionation via diffusion in the lower pH model runs than higher pH runs.
The choice of kinetic rate parameters and the associated isotopic kinetic fractionation factors for hydration and hydroxylation, as well as the fractionation of oxygen isotopes between water and hydroxide has a demonstrable effect on model parameters.In Figure S7, we test the effects of these choices on model output.These parameters were taken from Guo (27) (Fig. S7a, and Figure 3 of main text), Watkins and Devriendt (26) (Fig. S7b), and Chen et al., (66) (Fig. S7c).
All models subsequently use the same intrinsic clumped isotope fractionation factors, which are taken from Guo (27), as implemented by Watkins and Devriendt (26).Changing these parameters appears to result in different degrees of departure from ∆47/∆48 equilibrium, however with very similar overall behaviors.In the attached code (Data S2), it is possible to change between these parameters in the loadEqs.mscript.

Testing for possible isobaric effects on ∆47/∆48 by sulfur and oxidizable contaminants:
Isobaric interference is a problem in natural samples with mixed compositions, as organic components and minerals can release gases during acidification that would interfere at the 44-49 amu range measured in CO2 for clumped isotopes.Sulfide minerals are a common accessory mineral in methane seep deposits, as the sulfide derived from SD-AOM interacts with metal ions (e.g., Fe) to form minerals such as pyrite.The oxidation of such minerals during laboratory handling and analysis could potentially release sulfur oxide, whose most common isotopologue is 32 S 16 O, with a mass of 48.Thus, we deliberately sought out to test if this has had any measurable effect on the clumped isotope composition of our samples.This was done with two experiments, 1) deliberate contamination of a pure carbonate sample, and 2) bleaching of methane seep carbonate samples.This second experiment also serves as a test for other sources of isobaric interference as well.
A sample of Carrara marble was deliberately contaminated with 2 wt.% pyrite to test for these effects.Comparison between the contaminated and uncontaminated samples (shown in table 1 of the main text) showed no significant difference in isotopic composition, indicating that the presence of pyrite does not produce a measurable effect on the clumped isotope composition.
It has also been postulated that the presence of organic matter in carbonate samples may evoke a bias in clumped isotope values if not being removed prior to phosphoric acid digestion (e.g., Bergmann et al., 2018).To test if oxidative cleaning of samples has an effect, aliquots of Formosa Ridge carbonate were bleached for 24 hours in a 3wt.% sodium hypochlorite solution.
∆47 and ∆48 values of treated and untreated samples, however, were indistinguishable with 2SE, demonstrating that these samples are unlikely to be affected by isobaric contamination (Table 1).

Expressing isotopologue models' output as dδ/dt or d∆/dt
The numerical models presented by Chen et al.( 66), Uchikawa et al. (67), Watkins and Devriendt (26), Guo (27), and others calculate the instantaneous rate of change of the many isotopologues of dissolved inorganic carbon.Performing these calculations while tracking isotopologues as independent species makes calculation of isotopologue mixing and kinetic isotope effects much simpler.However, when interpreting these results, it is often helpful to recalculate the output values in delta-notation.The math governing this has been presented before with the following equations (where the isotopologue weight of a carbonate ion is given in brackets), in this case omitting the contribution of 17 O isotopologues at each mass following the conventions of Chen et al.( 66), Uchikawa et al. (67) and Watkins and Devriendt (26).
− 1000 (Equation S3) If the rate of change for each isotopologue (written here as [X]') can be calculated using a model, for instance as shown for mass-dependent diffusive fractionation earlier in this supplementary file, then the rate of change of δ 13 C, δ 18 O, ∆63 and ∆64 (written here as δ 13 C' δ 18 O' ∆63' and ∆64') can be calculated using the Constant Rule, Constant Multiple Rule, Product Rule and Quotient Rule to the following derivatives.
Converting the ∆63 and ∆64 values from equations S5, S6 into the more conventionally displayed CO2 ∆47 and ∆48 values for calcium carbonate is performed with the empirical linear transformations from Fiebig et al., (34).
(Equation S13) The relative rates of change of these parameters give the slopes of the lines for diffusive fractionation, CO2 dehydration and CO2 dehydroxylation shown in figure 3a and 3b in the main text.
Isotopologue end-member mixing: what (equilibrium) end-members could account for the observed variation in δ 13 C, δ 18 O, ∆47 and ∆48?Non-linear mixing effects in clumped isotopes are predictable and mathematically well-described within the published literature having been explicitly described for ∆47 by Defliese and Lohmann (68).Several studies have since presented theoretical and experimental tests of mixing for the dual-clumped isotope system (33,38,69).These non-linear mixing effects are a consequence of how the ∆47 and ∆48 parameters are calculated, as they are defined relative to the stochastic abundance of isotopes, which is the product of two values, represented using the reduced isotopologue approximation used by Watkins and Devriendt this is shown as ∆47 = 1000 R47/(R45R46) -1000.Because mixing lines between two end-members with different R45 and R46 values will have calculated stochastic R47 values that are parabolic in shape, whereas the mixing line of measured R47 values will be linear, there is an offset in these lines that appears to be a non-linear effect when only the ∆47 value is plotted, even though each individual isotopologue as behaves as a conservative "linear" mixing system.Whereas ∆47 mixing effects can be either positive or negative, ∆48 only ever exhibits positive mixing effects.It is possible to tune endmember δ 13 C and δ 18 O values to produce any measured slope of ∆47 and ∆48 values, such an exercise can also yield a minimum degree of internal heterogeneity necessary to produce a measured effect (38,69).
The variance in δ 13 C and δ 18 O values of carbonates at Site F have given previous workers reason to expect variable sources of carbon and water.Determining the possible contribution of endmember mixing at Site F provides a simple hypothetical question: can end-member mixing produce the observed disequilibrium effects in carbonate ∆47 and ∆48 values at Site F? In this exercise, it is necessary to assume that ∆47 and ∆48 values are at equilibrium, as not assuming this would necessitate the invocation of further kinetic isotope effects, which would negate the hypothesis.We further assume the equilibrium (clumped) temperatures of both endmembers are identical, that is to say heterogeneity is driven by the bulk composition of the end-members, this is done to minimize the degrees of freedom in our model.Because the observed ∆47 values are more negative, the δ 18 O and δ 13 C values must be anti-correlated with one another (68).To yield the observed slope of ∆47 and ∆48 values, the difference in end-member δ 18 O values must be roughly -0.2 to -0.6 the difference in end-member δ 13 C values.From this, we modeled mixtures with ∆δ 13 C values of -25‰ from a starting δ 13 C value of -45‰ and differences in δ 18 O of -5‰ -10‰ and -5‰ from an equilibrium value of 3.5‰.
In Figure S8, we show end-member mixing models that would be sufficient to describe the measured seep composition.The δ 13 C values of the modeled end-members are reasonable for methane-derived carbonate in this region, which has been measured at -70‰ (53).These mixing lines are plotted relative to the measured data in Figure S8.These models produced a sufficient mixing effect to describe the observed dataset, and necessitate an end-member water composition approximately 10‰ enriched relative to seawater.Experimental work at 4°C has shown that water in methane hydrate is enriched by 1.7‰-2.7‰relative to surrounding water (19), thus methane hydrate-associated water is not a plausible source of such profoundly δ 18 Oenriched water.Clay mineral dehydration and mixing with meteoric water both constitute endmembers with sufficiently dissimilar δ 18 O to seawater to account for such mixing effects (70); we conclude that these mixing effects are a factor worth considering, particularly in systems where δ 18 O can be extremely variable, such as with freshwater-seawater mixing zones.
In the analysis of the resulting carbonate systems, these mixing effects could be resolved in several circumstances, A) if carbonate derived from both end-members is present in the sediment, which is then homogenized and analyzed.B) if DIC from both end-members is mixed and subsequently mineralized more rapidly than the time required for the isotopic reequilibration of the DIC.(66), and the associated references therein.Analogue for seawater and freshwater mixing zones with methane-derived carbon more abundant in the freshwater endmember.
-70 -60 -50 -40 -h of the main text.By varying the methane flux, the exact magnitude of this offset can vary, as shown in Figure 3c-d of the main text.Testing the effects of other parameters on the COAD-MS model: The models discussed in the main text were implemented with model conditions that represented the present day at Site F, (i.e.Temperature = 3.5°C, salinity = 35 ppt, [Ca 2+ ] = 10 mM, pH = 7.7).Testing the model's behavior at other temperatures, salinities, ion concentrations and pH values would give predictions for how these systems would behave in other locations and seawater chemistries.Here, we re-run the model used to generate Figures 3g and 3h of the main text, with different parameters.

Fig. S1 .
Fig. S1.Alkalinity/DIC/pH/Wcalcite plot with vectors showing the net effect of aerobic oxidation of methane (AeOM), anaerobic oxidation of methane (AOM), aerobic oxidation of sulfide, and carbonate precipitation.Color map shows pH of solution as calculated using CO2SYS(57).Dashed contours show saturation state of calcite, as calculated using CO2SYS(57), assuming unchanging seawater Ca 2+ concentration of 10mM.All calculations use a salinity of 35ppt at 3.5°C at 1100 meters below sea surface.
Model runs identical to Figure3g,h of the main text, but with varying model temperature.
Model runs identical to Figure3g,h of the main text, but with variable seawater boundary condition pH.
Model runs identical to Figure3g,h of the main text, but with different values for k+1, k+4 (and associated kinetic isotope effects), as well as the aOH-H2O.Parameters taken from Guo(27),Watkins and Devriendt(26),Chen et al.,

Fig. S8 .
Fig. S8.Three mixing lines between the extrapolated equilibrium end-member, and hypothetical end-members with anti-correlated differences in δ 13 C and δ18 O, but identical ∆47 and ∆48 values.

Fig. S9 .
Fig. S9.Three mixing models between endmembers with correlated δ 13 C and δ 18 O values.Analogue for seawater and freshwater mixing zones with methane-derived carbon more abundant in the freshwater endmember.

Table S1 .
Default composition of initial model state, and boundary condition seawater flux during model run.Excel table containing raw data for clumped isotope analyses used in this study, alongside anonymized unrelated data that was used for error minimization following Daëron (58) and Bernecker et al. (55)..zipfile containing codes, data, and models used to generate Figure 2 and 3 of the main text.Also available via the following Zenodo DOI.https://doi.org/10.5281/zenodo.10954482